rm(list = ls(all = TRUE)) #clear workspace

setwd("~/Dropbox/REDISTRICTING/SIMULA_MVN")
source("02_get.parms.R")
source("03_simulate.R")
set.seed(12345) #Random seed for replication

# Computos 2012
data <- read.table("DIPUTADOS_2012.txt", header=TRUE, sep="\t")
coal <- data$COALICION.PRI.PVEM != 0
no.coal <- data$COALICION.PRI.PVEM == 0
data <- with(data, 
             data.frame(
               PAN=PAN,
               PRI=PRI + coal*(PVEM + COALICION.PRI.PVEM),
               Mov.Prog=PRD+PT+MOVIMIENTO.CIUDADANO+
                 COALICION.PRD.PT.MC+COALICION.PRD.PT+COALICION.PRD.MC+COALICION.PT.MC,
               OTH=NUEVA.ALIANZA + (PVEM*no.coal)))

# Curules observadas
obs.seats <- apply(data, 1,function(x) max(x)==x)
obs.seats <- rowSums(obs.seats)

# Simula escenario observado
#----------------------------
# Parametros
parms <- get.parms(data)

# Simula
sim.obs <- replicate(1000, simulate(n=300, mu=parms$mu, Sigma=parms$sigma), simplify=FALSE)

sim.out <- do.call(rbind, sapply(sim.obs, "[", "out"))
sim.vote <- do.call(rbind, sapply(sim.obs, "[", "mean.vote"))
sim.seats <- do.call(rbind, sapply(sim.obs, "[", "seats"))

# Compara con resultado observado
rbind(obs=obs.seats, apply(sim.seats, 2, quantile, c(0.5,0.025,0.975)))

# Simula escenario hipotetico
#----------------------------
hip.mu <- c(.31,.31,.31,.07)
hip.mu <- log(hip.mu[-length(hip.mu)]/hip.mu[length(hip.mu)])
#hip.mu <- parms$mu
#hip.mu[1:3] <- mean(hip.mu[1:3])

#Variar mu, misma sigma ("swing ratio")
sim.hip <- replicate(1000, simulate(n=300, mu=hip.mu, Sigma=parms$sigma), simplify=FALSE)

sim.hip.out <- do.call(rbind, sapply(sim.hip, "[", "out"))
sim.hip.vote <- do.call(rbind, sapply(sim.hip, "[", "mean.vote"))
sim.hip.seats <-do.call(rbind, sapply(sim.hip, "[", "seats"))

# Curules
apply(sim.hip.seats, 2, quantile, c(0.5,0.025,0.975))

##Ternary plots
##-------------
library(gridExtra)
library(ggtern)
# Density
plot.dens <- stat_density2d(method = "lm", fullrange = T, 
                            n = 200, geom = "polygon", 
                            aes(fill = ..level.., 
                                alpha = ..level..))
# Win-set
lines <- data.frame(x = c(1, 0, 1)/3, 
                    y = c(1, 1, 0)/3, 
                    z = c(0, 1, 1)/3,
                    xend = c(1, 1, 1)/3, 
                    yend = c(1, 1, 1)/3, 
                    zend = c(1, 1, 1)/3)
win.set <- geom_segment(data = lines, 
                        aes(x, y, z, 
                            xend = xend, yend = yend, zend = zend), 
                        color = "gray48", size = 0.7)

## Plot - Resultados distritales
png("plot1_observado.png")
  ggtern(data = data[,1:3], aes(PAN, PRI, Mov.Prog)) +
  theme_bw() +
  geom_point(color="blue", alpha = 0.5) + 
  geom_confidence(color="red") +
  labs(title = "Voto distrital observado") + 
  guides(color = "none", fill = "none", alpha = "none") +
  win.set
dev.off()

# Plot - Simulated election
png("plot2_simulado.png")
  ggtern(data = as.data.frame(sim.out[,1:3]), aes(PAN, PRI, Mov.Prog)) + 
  theme_bw() +
  plot.dens + scale_fill_gradient(low = "blue",high = "red") +
  labs(title = "Voto distrital simulado") + 
  guides(color = "none", fill = "none", alpha = "none") +
  win.set
dev.off()

# Plot - Escenario hipotético
png("plot3_hipotetico.png")
  ggtern(data = as.data.frame(sim.hip.out[,1:3]), aes(PAN, PRI, Mov.Prog)) +
  theme_bw() +
  plot.dens + scale_fill_gradient(low = "blue",high = "red") +
  labs(title = "Escenario hipotético") + 
  guides(color = "none", fill = "none", alpha = "none") +
  win.set
dev.off()

